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Neural mass signals from in-vivo recordings often show oscillations with frequencies 
ranging from <1 to 100 Hz. Fast rhythmic activity in the beta and gamma range can be 
generated by network-based mechanisms such as recurrent synaptic excitation-inhibition 
loops. Slower oscillations might instead depend on neuronal adaptation currents whose 
timescales range from tens of milliseconds to seconds. Here we investigate how 
the dynamics of such adaptation currents contribute to spike rate oscillations and 
resonance properties in recurrent networks of excitatory and inhibitory neurons. Based 
on a network of sparsely coupled spiking model neurons with two types of adaptation 
current and conductance-based synapses with heterogeneous strengths and delays we 
use a mean-field approach to analyze oscillatory network activity. For constant external 
input, we find that spike-triggered adaptation currents provide a mechanism to generate 
slow oscillations over a wide range of adaptation timescales as long as recurrent 
synaptic excitation is sufficiently strong. Faster rhythms occur when recurrent inhibition is 
slower than excitation and oscillation frequency increases with the strength of inhibition. 
Adaptation facilitates such network-based oscillations for fast synaptic inhibition and leads 
to decreased frequencies. For oscillatory external input, adaptation currents amplify a 
narrow band of frequencies and cause phase advances for low frequencies in addition 
to phase delays at higher frequencies. Our results therefore identify the different key roles 
of neuronal adaptation dynamics for rhythmogenesis and selective signal propagation in 
recurrent networks. 



Keywords: spike frequency adaptation, adaptation, oscillations, rate models, network dynamics, Fokker-Planck, 
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INTRODUCTION 

A prominent characteristic of cortical activity is its rhythmicity 
as shown by electroencephalography or the local field potential. 
Dominant oscillation frequencies in these signals range from < 1 
to 100 Hz and reflect synchronous activity of populations of neu- 
rons. Such oscillations are linked to behavioral states (Wang, 
2010) and involved in a variety of cognitive functions (Engel et al, 
2001; Fries, 2001; Melloni et al, 2007; Ghazanfar et al, 2008; 
Wang, 2010) as well as pathological conditions (Hammond et al, 
2007; Zijlmans et al, 2009; Uhlhaas and Singer, 2010). It is there- 
fore important to understand the mechanisms of oscillations in 
neuronal networks, how they are initiated and terminated, and 
how their frequency is determined. 

Fast rhythmic activity in the beta and gamma band (>20 Hz) 
can be generated by network-based mechanisms, such as synap- 
tic excitation-inhibition loops or by feedback inhibition alone 
(Isaacson and Scanziani, 2011). In these scenarios the oscilla- 
tion frequency is largely determined by the inhibitory decay 
time constant (Brunei and Wang, 2003; Tiesinga and Sejnowski, 
2009). Low-frequency oscillations, on the other hand, could 
depend on slow transmembrane outward currents (Compte et al, 
2003; Gigante et al, 2007b; Destexhe, 2009), which are medi- 
ated by low-threshold voltage-dependent muscarinic (M) and 
high-threshold calcium-gated afterhyperpolarization (AHP) K + 



channels, respectively (Brown and Adams, 1980; Connors et al., 
1982; Stocker, 2004). These currents cause spike frequency adap- 
tation and are typically more pronounced in cortical regular 
spiking pyramidal (excitatory) neurons compared to fast spiking 
(inhibitory) interneurons (La Camera et al., 2006). Both, the M 
and AHP type K + currents, are susceptible to cholinergic mod- 
ulation (McCormick, 1992). Their kinetic time constants range 
from milliseconds to seconds (Abel et al, 2004; Manuel et al, 
2005) and can be pharmacologically manipulated (Pedarzani 
etal., 2001). 

Here we study the interplay of the dynamics of such adapta- 
tion currents with synaptic excitation and inhibition in recurrent 
networks of excitatory and inhibitory neurons. Specifically, we 
ask (1) how adaptation can generate slow oscillations, (2) how it 
modulates faster rhythms based on synaptic interaction, and (3) 
how adaptation affects resonance properties of the network. 

In-vivo recordings from behaving animals have revealed that 
even when the population activity oscillates, the spike trains of 
the constituent neurons are rather irregular and display Poisson- 
like characteristics (Fries, 2001; Wang, 2010). This stochasticity 
in neuronal responses allows us to derive a mean-field model 
from a recurrent network of adaptive spiking model neurons 
coupled through conductance-based synapses with heteroge- 
neous strengths and delays. Our approach is based on the 
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Fokker-Planck (FP) formalism (Brunei, 2000; Deco et al., 2008) 
and efficiently describes the activity of large networks where the 
features of the spiking neurons (i.e., the model parameters) are 
retained. Using this method we analyze network responses to 
constant as well as rhythmic external input. In particular we 
describe asynchronous irregular states with constant steady-state 
activity as well as oscillatory states and their properties. We val- 
idate our mean-field results qualitatively by large-scale network 
simulations. 

METHODS 

We first describe our network model containing two popula- 
tions (excitatory and inhibitory) of adaptive spiking neurons 
with delayed conductance-based synaptic coupling. Based on that 
model we then derive mean-field model equations and solve them 
numerically to obtain distributions of the membrane potentials 
and instantaneous spike rates. 

NETWORK MODEL 

We consider a network of N = Ns + Nx adaptive exponen- 
tial integrate-and-fire neurons (aEIF) proposed by Brette and 
Gerstner (2005), where Ng and Nj are the numbers of excita- 
tory and inhibitory neurons, respectively. The dynamics of the 
i-th neuron of population a e {£, X] is described by 

C^=WV?)-<+^ y n,(^0 (1) 

with reset condition 

f V a ■= Vr 

ifV? > Vcutthen ' T (3) 
' \wf := + b. 

The first Equation (1) is for the membrane potential V a , where 
the capacitive current through the membrane with capacitance C 
equals the sum of ionic currents Ii 0D , the adaptation current wf 
and the synaptic current i^L ; . The ionic currents are given by 

v-v T 

I io n(V):=g L (E L -V)+g L A T e a t , (4) 

where the first term on the right-hand side describes an Ohmic 
leak current with conductance gi and reversal potential Ei. The 
exponential term with threshold slope factor At and threshold 
potential Vj approximates the Na + -current which is responsi- 
ble for the generation of spikes, assuming that the activation of 
Na + -channels is instantaneous and neglecting their inactivation 
(Fourcaud-Trocme et al., 2003). Equation (2) governs the dynam- 
ics of the adaptation current w a , where r w denotes the adaptation 
time constant and a quantifies a conductance that mediates sub- 
threshold adaptation. A spike is said to occur at the time when 
V" diverges to infinity, but in practice a finite "cutoff" value 
Vcut is chosen. When V a crosses V cut from below, V a is set to 
the reset potential V, and w a is incremented by b, cf. condi- 
tion (3). In this way spike-triggered adaptation is included in the 



model. Immediately after the reset, V a and w a are clamped for a 
refractory period T re f. 

The aEIF model has been shown to reproduce a broad range of 
subthreshold dynamics (Touboul and Brette, 2008) and spike pat- 
terns of cortical neurons (Naud et al., 2008) and can well predict 
their spike times (Jolivet et al., 2008) and post-stimulus time his- 
tograms (Pospischil et al, 2011). Importantly, the subthreshold 
and spike-triggered adaptation components of this model have 
been shown to capture the effects of the M and AHP currents 
in a detailed biophysical neuron model, respectively (Ladenbauer 
et al., 2012). 

Neuron i of population a receives total synaptic current 

; ; ; 

which is the superposition of synaptic inputs J"' 6 from K ext 
external excitatory neurons, from Kg excitatory neurons of 
the network and from Kx inhibitory neurons of the network. 
j is the index of the respective presynaptic neuron. The synaptic 
current /™' Y caused by neuron j of population y e {ext, £, X} is 
modeled using delta functions, 

q<*<yf, t) := c;«- ext J2 & ( f - ff) ( E £ ~ K) ( 6 ) 

k 

£P(V?, t) := Cjf £ 5(f - tf - 4 P ) (Ep - V?) , (7) 

k 

where |3 e {£, 1} denotes the presynaptic population. /*' Y are 
dimensionless synaptic efficacies drawn from a Gaussian distri- 
bution with mean J a y and standard deviation AJ a y . Here we 
consider that / a Y = J Y and AJ a y = A/ y depend only on the 
presynaptic population y. f| is the /c-th spike time of neuron j 
from the respective population. Eg and Ex denote the excitatory 
and inhibitory reversal potentials, respectively, cf*^ is the synaptic 
delay, sampled using a bi-exponential probability density 

1 / d-dp d — dq \ 

P?W.= — [e- * ~e * j (8) 

for positive delays d, where c?o is the minimal delay and zy, axe 
the rise and decay time constants, for each pair of populations. 
In the model we use two different delay distributions p £ d and 
which do not depend on the postsynaptic population as for the 
synaptic weights. For a schematic diagram of the network, see 
Figure 1 . 

We assume the neurons from the external population gener- 
ate spike times according to Poisson processes with rates r" xt (f). 
The spike rate of each population a e {£, X) at time t is given by 
the average number of spikes of neurons from the corresponding 
population in the interval [f, t + At], 

/= 1 k 
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FIGURE 1 | Network architecture. Each of N £ excitatory and N x inhibitory 
neurons receives excitatory input from Kext external neurons with mean 
synaptic strength J ext as well as synaptic input from Kg (K x ) excitatory 
(inhibitory) neurons of the network with mean strength J e (J x ) and delays 
distributed according to p £ d (pj). 



In the mean-field limit N — > oo, At — > 0 we obtain a continuous 
population spike rate r a (t) (see below). 

We selected the following parameters for the neuron model: 
C = 200 pF, g L = 10 nS, E L = -70 mV, A r = 1 mV, V T = 
-50 mV, V r = -70mV, V cut = -40mV, and T re f = 1.4 ms 
(Badel et al, 2008; Destexhe, 2009). For excitatory neurons the 
adaptation parameters were varied within reasonable ranges: 
r vv € [5, 1000] ms, a e [0, 10] nS, b e [0, 50] pA. For inhibitory 
neurons adaptation was neglected (a = b = 0) since it was found 
to be weak in fast spiking interneurons compared to pyramidal 
cells (La Camera et al., 2006). 

The network parameter values were N £ =40,000, Nj = 
10,000, K ext = 1600, K £ = 1600, K T = 400, E £ = OmV, E T = 
-80 mV, / ext = 0.003, J £ = 0.003, and A/ y = 0.1/ Y with y e 
{ext, £ , T] (Brunei and Wang, 2003). To adjust the balance of 
recurrent synaptic excitation and inhibition we introduce the 
parameter 



h\E L -E T \ 
J £ \El-E £ \' 



such that r £ = rj in case of uncoupled populations of neurons, 
i.e.,/ £ =J X = 0. 

MEAN-FIELD MODEL 

We reduce the two-population network of aEIF neurons to the 
mean-field model in three steps. First, we replace the synap- 
tic current fluctuations by a Gaussian white noise process via 
the diffusion approximation. Next, we take a mean-field limit 
to formulate the stochastic network model in terms of two cou- 
pled deterministic scalar partial differential equations (PDE). 
Finally, to allow for efficient numerical computation we reduce 
the number of variables in these equations using an adiabatic 
approximation. 

Diffusion approximation 

We approximate the total synaptic current 7^ n ; of Equation (5) 
by its mean plus a fluctuating Gaussian part, which is jus- 
tified by the following physiologically plausible assumptions: 
(1) The number of synaptic inputs to a neuron is large, i.e., 
K eil t, K £ , Kj ^> 1 (Destexhe et al., 2003) and (2) the postsynaptic 
potential amplitudes elicited by individual presynaptic spikes are 
smaU, i.e., J m \E £ - V\, } £ \E £ - V\, h\Ej - V\ « V cut - V r 
(Williams and Stuart, 2002). We further assume that (3) the 
network connectivity is random and sparse, i.e., K £ , Kj <^ N, 
and that (4) presynaptic spike times are represented by Poisson 
processes which are homogeneous in each small time interval. 
The total synaptic current can then be written as (Brunei, 2000; 
Nykamp and Tranchina, 2000; Renart et al., 2004; Richardson, 
2004; Gigante et al, 2007b) 



Cn,i ~ ^a,i(Vf , f) + a aj! (Vf , t)r\ { (t), 



(11) 



(10) 



where [i a i and cr a , are the infinitesimal mean and standard 
deviation of 7^ n ; , respectively, and r|, is a Gaussian white 
noise process with 8-autocorrelation. The infinitesimal mean is 
given by 



which is the ratio of total charges induced at rest (Kumar et al., 
2008). g determines Jx and thus A/x for fixed } £ and was var- 
ied in [0.8, 4] which yields a physiological range of inhibitory 
postsynaptic potential amplitudes (Tamas et al, 1997). Note 
that the value of g that corresponds to balanced mean recur- 
rent excitatory and inhibitory synaptic currents depends on the 
mean membrane potential for each population. The effect of 
a spike of presynaptic neuron j on neuron i is mediated by 
a delayed instantaneous increment or decrement of the post- 
synaptic membrane potential, cf. Equations (1), (5), and (7). 
This implies that reflects the conduction delay as well as 
delays in the synaptic kinetics. We therefore chose the parame- 
ter values of p £ d and such that conduction delays as well as 
typical time courses of excitatory AMPA and inhibitory GABAa 
synaptic receptors are taken into account. The values we selected 
were d 0 = 1ms, xf e [1.25, 1.5] ms, tf e [1.5, 2]ms, r, x e 
[0.55, 1.25] ms, and rj e [1.5, 5] ms. The input rate of the exci- 
tatory population rf^ was varied in [1, 12.5] Hz. r^ xt was chosen 



with 



li aA := lim J ^ L (12) 



At 



^ai = C(J5p-V?)/pKj,(r|j*pp)(0, 



(13) 
(14) 



where (•) denotes the expectation operator. The infinitesimal 
variance is 



a : = lim 



(// +At W^) 2 ) + °( Af2 ) 



Al 



(15) 



(0 2 + (<<) 2 + (<,) 2 
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with 



°Z = C( - E £ ~ V DV + A/it) WextC) 



C(£p - Vf)J(jj + A/ 2 ) %(r p * p p )(0, 



(16) 



(17) 



where P € {£ ,T) and * denotes convolution. In Equations (13), 
(14), (16), and (17) we have used that the presynap- 
tic Poisson processes, the synaptic weights and delays are 
independent. 

Mean-field limit 

We analyze networks of sparsely coupled neurons, i.e., the prob- 
ability for a connection between any pair of neurons is low, 
cf. assumption (3) above. For large N correlations between the 
fluctuations of synaptic currents of different neurons become 
negligible, i.e., (r\i(t)r\j(t)) = 0 for i ^ j. In the mean-field limit 
N — » oo the network model Equations (1)— (4), Equations (11)- 
(17) can be described by two FP equations — one for each popu- 
lation a — which are delay-coupled by the population spike rates 
r£ and rx, 



dt 



+ 



9^ 
dV 



+ 



dS™ 
dw 



with 



kon{V) ~ W+ [Xq 

C 

2C 2 dV 
a{V — Ei) — w 

Pa- 



0 



a a da a 



(18) 



-^—)p a (19) 
2C 2 dV ' 



(20) 



p a (V, W, t) is the probability density to find a neuron of popula- 
tion a in the state (V, w) at time t. S^(V, w, t) and S™(V, w, t) are 
the probability fluxes in positive V and w direction, respectively. 
Note that we used the Stratonovich interpretation of the under- 
lying stochastic equations (Risken, 1996; Richardson, 2004). To 
account for the reset condition (3) the flux through the cutoff 
voltage V cut at w is re-injected after the refractory period T re f at 
V r , w + b, i.e., 



The spike rate of population a is given by the integral of the cutoff 
fluxes, 



r a (t) 



f s v a 
Jr 



(Vcut, w, t)dw. 



(25) 



At any timepoint t the histogram of the membrane potentials of 
neurons in population a can be seen as a sample drawn from 
the probability density p a (V, t) which is governed by the FP 
equation. 

Adiabatic approximation 

Solving the 2+1 dimensional PDE (Equations 18-20) with cor- 
responding reset and boundary conditions (2 1 )-(24) numerically 
is possible but computationally demanding. We therefore reduce 
the dimensionality of the FP system Equations (18)-(20) assum- 
ing the timescales of membrane voltage and adaptation current 
dynamics are separable. This is justified by the observation that 
the dynamics of neuronal adaptation is significantly slower than 
the other in the model system such as membrane time con- 
stant and average inter-spike interval (Womble and Moises, 1992; 
Stocker, 2004). Under this assumption, the adaptation current 
of each neuron can be seen as an efficient integrator that filters 
the fluctuations in the neuronal activity. We approximate w a (t) 
in Equation (2) by its population average w a (t), which evolves 
according to 



dw a 

Xv, -fc - a (( v )pa(V,t) 



E L ) - w a + r w br a (t), 



(26) 



where (-) p denotes the average over the density p (Brunei et al., 
2003; Gigante etal, 2007b). The probability density p a (V, t) then 
satisfies the 1 + 1 dimensional FP equation 



V 



dpa dS± 
dt dV 



(27) 



where again S„ is the probability flux defined in Equation (19) 
and w := w a (t) appears as a system parameter. The reset condi- 
tion is 



lim S*(V, t) - lim SUV, t) 

V\V r V\V r 



(28) 



T ref ). 



lim SUV, w+b,t)- lim SUV, w + b,t) 

V\V r VfV r 

= S^(V cut , w, t - T ref ) Vw e R. 



(2\) and the boundary conditions (23)-(24) become 



This implies that in general p a is not differentiable at the 
line V = V r . The boundary conditions are reflecting for 
w —>■ ±oo, V —>■ —oo and absorbing for V = V cut , 



lim S£(V, w) = 0 We (-oo, Vc Ut ] 
lim SY(V, w) = 0 W e R 

Pa(V cut , W) = 0 Vw € R 



lim SUV) = 0, 

V^— oo 

Pa (V cut ) = 0. 



(29) 
(30) 



The population spike rates are given by the corresponding fluxes 
through the cutoff voltage, 



r a (t) = SUV cut ,t). 



(31) 



(22) 
(23) 

Note that the adiabatic approximation described above could be 
(24) applied repeatedly for additional slow variables. 
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NUMERICAL SOLUTION 

We solved the reduced FP Equation (27) subject to con- 
ditions (28)-(30) and mean adaptation current dynamics 
(Equation 26) forward in time until either steady states r|? 
withr^° := lim^oo r a (t) or stable oscillatory states were reached. 
The probability densities pg, px were initialized using normal- 
ized Gaussians with mean 0.5 • (V r + Vj) and standard devi- 
ation 0.2 • (Vj — V r ). We applied a first-order finite volume 
method on a finite and non-uniform grid Vq < V\ < ■ ■ ■ < Vjv v 
using upwind-fluxes to stabilize the numerical solution (LeVeque, 
2002). Time was discretized using the implicit Euler method 
on an equidistant grid, i.e., f„ + i — t n = At. The resulting lin- 
ear equation systems were solved with a preconditioned Krylov 
subspace method in each time step. Specifically, BiCGSTAB 
(van der Vorst, 1992) was used in combination with an incom- 
plete LU decomposition preconditioner (Saad, 2003) that strongly 
improved the convergence speed. 

W£ was initialized with values w^(0) 6 [0, 500] pA (and 
wj = 0). The other parameters were Af = 50 |xs, min m A V m = 
1 (jlV with AV m := V m+1 - V m , V 0 := -100 (jlV, V Nv = V cut 
and N v = 256. 

We complemented the mean-field results with numerical sim- 
ulations of the network model Equations (l)-(4) using a Runge- 
Kutta second order method implemented in Brian 1.4 (Goodman 
and Brette, 2009) with a time step of 50 \is. 

In case of stable periodic population spike rates the oscilla- 
tion frequency was determined by the dominant frequency of the 
Fourier spectrum of rg over the last 2 s of runtime. 

RESULTS 

ADAPTATION MEDIATES OSCILLATIONS 

To examine how the interplay of adaptation and recurrent synap- 
tic input shapes network dynamics we vary the type, strength and 



timescale (parameters a, b, and r w ) of adaptation for excitatory 
neurons as well as the strength of synaptic inhibition (parame- 
ter g) across networks. Adaptation currents are disregarded for 
inhibitory neurons, which is supported by experimental obser- 
vations, see the section Methods. We consider constant rates 
r fxt> 'rat f° r me external Poisson-inputs and identical delay dis- 
tributions p A = pj. First, we examine steady-state spike rates, 
oscillation amplitudes and frequencies for networks with differ- 
ent values of spike-triggered adaptation b and inhibition strength 
g, see Figure 2A. All networks without adaptation (a = b = 0) 
settle into asynchronous states with constant population rates 
that decrease with increasing g. For networks with increased b 
slow oscillatory states become stable if recurrent excitation is suf- 
ficiently strong. The larger b is, the less recurrent excitation is 
necessary for sustained oscillations. Amplitude and period of the 
oscillatory rate decrease with an increase of b and g, respectively. 
Thus, in networks where recurrent synaptic excitation dominates 
inhibition at least slightly, spike-triggered adaptation b generates 
spike rate oscillations. The dynamics of an example network is 
shown in Figure 2B. The evolution of the population spike rates 
rg,rx, membrane potential probability densities pg , px and adap- 
tation current W£ display periodic bursts of population activity. 
As a validation of the findings above using the mean-field model 
the activity of simulated large networks of spiking neurons is 
shown in Figure 2C. The raster plots reveal population bursts 
when b is increased and g is small. An asynchronous state with 
low population activity occurs if g is increased. If in addition 
adaptation is removed (a = b = 0) the network settles into an 
asynchronous state with increased spike rates. 

The mechanism that generates these oscillations is a loop of 
recurrent excitation, build up and decay of adaptation current as 
indicated in Figure 2B. A low level of population activity is ini- 
tiated by the external input rf xt and recurrent synaptic excitation 




FIGURE 2 | Population bursts caused by spike-triggered adaptation. 

(A) Top: Spike rate r £ of the excitatory population as a function of the 
strength of inhibition g for networks without spike-triggered adaptation 
(b = 0, black) and with increased levels of b (0.025 nA, brown and 0.05 nA, 
red). In case of stable oscillatory states the maxima and minima of the 
periodic re are shown by dashed lines. Solid lines represent asynchronous 
states. Arrows indicate balance of recurrent excitation and inhibition for 
both populations. Bottom: Corresponding oscillation frequencies 1. 
t w = 200 ms, a= 0, and rf xt = 6.25 Hz. The parameter values for both delay 
distributions p e d , were %, = 1.5ms and x t j = 2ms. For other model 



parameters see the section Methods. (B), Top: Time-dependent spike rates 
ffiCO (green) and r x (f) (orange, dashed) for the parameter values 
6 = 0.05nA and g='\, as indicated in (A) by red dots. Center: 
Corresponding membrane potential density p £ (V, f). Bottom: Corresponding 
mean adaptation current wj(f). (C): Raster plots of simulated networks of 
N = 50,000 aEIF neurons for b = 0.05 nA, g = 0.85 (top), b = 0.05 nA, 
g = 1.05 (center) and b = 0, g = 1 (bottom). The spike times of 200 
excitatory neurons and 50 inhibitory neurons, all randomly selected, are 
shown by green and orange dots, respectively. x„ = 200 ms, a = 0, and 
ff xt = 3.75 Hz. Other parameter values as in (A). 
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boosts the activity, thereby increasing the adaptation current W£ 
through b in a spike rate dependent way. The adaptation current 
in turn acts as a negative feedback which eventually outweighs 
the recurrent excitation. The population activity drops rapidly 
and the adaptation current decays slowly. Upon recovery from the 
adaptation current the cycle starts again. 

Next, we investigate how these oscillations are affected by 
the external input rf xt , the subthreshold adaptation conduc- 
tance a and the adaptation timescale r w , see Figure 3. The 
existence of adaptation-induced oscillations is quite sensitive to 
the level of rf^ (Figure 3A). Periodic activity is stable for small 
values of rf^ (above threshold). While oscillation frequencies 
increase monotonically with increasing rf^, oscillation ampli- 
tudes increase initially for a small interval of rf xt values and 
decrease over the following interval. For larger values of r^ xt 
oscillatory activity is destabilized and asynchronous states occur. 
Interestingly, an increase in a does not lead to oscillations. On 
the contrary, periodic population bursts are destabilized by a. 
The dependence of oscillation amplitude and frequency on r w 
is shown in Figure 3B. Stable oscillations exist for a large range 
of values of r w , where the frequencies decrease with increasing 
t w . Oscillations are unstable for small adaptation timescales in 
the range of the membrane time constant and for very large 
values of r w . 



ADAPTATION MODULATES FREQUENCIES OF NETWORK-BASED 
OSCILLATIONS 

Here we study the influence of adaptation on oscillations gen- 
erated by recurrent synaptic excitation-inhibition {£-!) loops. 
The pace of such oscillations is believed to be largely deter- 
mined by the decay of inhibition. To describe their dependence 
on the timescale of inhibition for various recurrent network 
regimes (from excitation dominated to inhibition dominated) we 
first consider networks of neurons without an adaptation cur- 
rent (a = b = 0), see Figures 4A,B- By varying the decay rj of 
inhibition and its strength (by parameter g) across networks we 
find that stable oscillatory states occur if inhibition is sufficiently 
slow in comparison to excitation. The oscillation frequencies 
increase with increasing external input spike rate rf^, increas- 
ing g and decreasing tJ, respectively. A low value of leads 
to frequencies in the low beta band (Figure 4A), for a higher 
value of rf xt the frequencies span the beta and low gamma bands 
(Figure 4B). Note that the network parameters can be adjusted to 
obtain higher oscillation frequencies. The generating mechanism 
underlying the oscillations is a loop of recurrent synaptic exci- 
tation and inhibition, initiated by the excitatory external input. 
We verified this by removing the recurrent excitatory input to the 
inhibitory population, which lead to a destabilization of the oscil- 
lations. For larger values of g as the ones used in Figure 4, the 




ext 



FIGURE 3 | Effects of subthreshold adaptation, external input, and 
adaptation timescale on population bursts. (A), Top: Spike rate r t - 
depending on the external input rf xt for networks without subthreshold 
adaptation (a=0, red) and with increased levels of a (5nS, violet and 
10 nS, dark blue). Maxima and minima of oscillating r £ are shown by 



dashed lines. Bottom: Corresponding frequencies f. fo=0.05nA, 
z w = 200 ms, g = 1, and other parameter values as in Figure 2A. (B): 
Maxima and minima of r e (top) and oscillation frequency as a function of 
the adaptation time constant x w . a = 0, /f xt = 6.25 Hz, and other 
parameter values as in (A). 
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A a. = 0 nS, 6 = 0 nA, 

r ext = 315 Hz 



dom. 



Inh. dominated 



B a = 0 nS, 6 = 0 nA, 
rf T , = 6.25 Hz 




ASYN 

IL^ 



a = 0 nS, 6 = 0.05 nA, D a = 5 nS, 6 = 0 nA, 




FIGURE 4 | Influence of synaptic inhibition and adaptation on 
network-based oscillations. (A-D): Existence of oscillatory states (OSC) 
and corresponding frequencies fas a function of the strength g and 
timescale of synaptic inhibition for networks with adaptation parameters 
and external input strengths as specified. Asynchronous states (ASYN) are 



indicated by white regions in the parameter space. Arrows mark balance of 
recurrent excitation and inhibition. On the left p £ d (green) and pj (orange) are 
shown for tJ = 1.5 ms, rj = 5 ms. xf was chosen such that the peaks of p^ 
and pj occur at the same delay value. t w = 200 ms, xf = 1.25 ms, and 
xS = 1.5 ms. For other parameter values see the Methods section. 



£-Z-loop mechanism is replaced by an Z-Z-loop that does not 
depend on recurrent excitation (not shown). Since adaptation is 
only exhibited by excitatory neurons, we disregard the parameter 
space where X-X-loop-based rhythmic activity occurs and focus 
on £ -Z-loop-based oscillations instead. 

An increase of spike-triggered adaptation or subthreshold cur- 
rent stabilizes oscillatory states also for faster recurrent inhibition, 
see Figures 4C,D. This change in single neuron dynamics causes 
oscillations in large parts of explored (g, rj)-space. In particular, 
for spike-triggered adaptation asynchronous states only occur in a 
small region of the parameter space. Interestingly, the oscillation 
frequencies are significantly reduced by either type of adaptation. 

Next, we investigate how the timescale of adaptation r w affects 
oscillations mediated by an £-Z-loop. In Figure 5A we show 
the dependence of amplitude and frequency of such oscillations 
on r w for networks with both adaptation components increased 
(a = 5 nS, b = 0.05 nA) and either dominant recurrent excitation 
(g = 1.05) or inhibition (g = 1.5). In both cases, stable oscilla- 
tory states exist for a large range of time constants. As r w increases 
the oscillation frequencies decrease while the amplitudes first 
increase abruptly and then decrease. The networks settle into 
asynchronous states for small r w (in the order of the membrane 
time constant) or large z w (several hundreds of milliseconds). 
Note that these effects of r w are similar if either a or b is increased 
individually (not shown). We validated these effects by simula- 
tions of aEIF neuron networks, see Figure 5B. The raster plots 
show that an increase in r w leads to a decrease in oscillation 
frequency and amplitude. 

ADAPTATION PROMOTES PERIODIC SIGNAL PROPAGATION 

To analyze how the resonance properties of recurrent networks 
in asynchronous states are influenced by adaptation currents, we 
here consider external Poisson-inputs with oscillatory rates with 
frequency/. Gain of input spike rate and phase difference between 
network and input spike rates as a function of input frequency 
for networks without (a = b = 0) and with adaptation (a = 5 nS, 



b = 0.05 nA) considering two adaptation time constants are pre- 
sented in Figures 6A,B. Excitation dominated networks without 
adaptation do not exhibit resonance at any frequency and show 
only phase delays. The presence of an adaptation current leads to a 
significant amplification of oscillations in the input which is par- 
ticularly strong at lower frequencies (of the beta band). This effect 
is pronounced for an increased adaptation timescale. In addition, 
adaptation causes a phase advance for low oscillation frequencies. 

In networks where recurrent inhibition dominates excitation 
on the other hand even in the absence of adaptation currents res- 
onance is shown for a high frequency band and phase advances 
for lower frequencies. Adaptation greatly enhances resonance and 
shifts the preferred frequency band to the high gamma range. 
The resonance effect is even stronger if the adaptation current 
is slower, i.e., r w increased. Although these effects of adaptation 
on resonance properties of recurrent networks are similar when 
either the subthreshold (a) or spike-triggered adaptation compo- 
nent (b) is increased individually, the dominant contribution to 
the frequency amplifications comes from b (not shown). We addi- 
tionally examined the response of single neurons to oscillatory 
noisy inputs using our mean-field model and found that adap- 
tation mediates resonance even in the absence of recurrent input 
(not shown). These results emphasize the importance of adap- 
tation for the amplification and thus propagation of oscillatory 
signals in neuronal networks. 

DISCUSSION 

In this work we have investigated the role of neuronal adaptation 
currents in shaping spike rate oscillations in large recurrent net- 
works of excitatory and inhibitory neurons. Based on a network 
of aEIF model neurons sparsely coupled through conductance- 
based synapses with heterogeneous delays and strengths driven by 
noisy external input, we used a mean-field method taking advan- 
tage of the FP equation. We simplified the problem by applying 
an adiabatic approximation and solved the resulting equations 
numerically. Using this method we obtain membrane potential 
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FIGURE 5 | Effects of adaptation timescale on network-based 
oscillations. (A) Top: Spike rate re as a function of the adaptation 
time constant r w for networks with dominant recurrent excitation 
(g=1.05, violet) and inhibition (a; =1.5, blue). Dashed lines indicate 
maxima and minima of oscillating r £ , solid lines represent constant 
f£. Bottom: Corresponding oscillation frequencies f. a=5nS and 



b=0.05nA. /f xt = 7.5Hz, zf = 1.25ms, r| = 1.55ms, = 0.98ms, 
and rj = 2ms. Other parameters as in Figure 4. (B): Raster plots of 
simulated networks of size N= 50,000 with g=1.5 and t w = 100ms 
(top) as well as t w = 400 ms (bottom), showing the spike times of 
200 excitatory and 50 inhibitory aEIF neurons. Other parameter values 
as in (A). 



distributions and population averages of spike rates and adapta- 
tion currents. At the same time, the dynamical properties of single 
neurons, i.e., the neuron model parameters, are retained in the 
derived mean-field network model. 

Alternative mean-field methods have been developed for 
conductance-based model neurons (Robinson et al, 2008) and 
recurrent networks thereof in asynchronous states (Shriki et al., 
2003), where spike rates are obtained without having to solve a 
PDE. Our approach based on the FP equation on the other hand 
treats noise in the synaptic inputs in more detail and allows for 
the calculation of membrane potential distributions in addition 
to spike rates. 

We chose the aEIF model because it provides a rich yet low- 
dimensional description of neuronal dynamics and includes a 
proper phenomenological description of the M and AHP adapta- 
tion currents. The effects of subthreshold (a) and spike-triggered 
adaptation (b) on response properties of aEIF neurons (measured 
by spike rate-input current relationships and phase response 
curves) match those of M and AHP adaptation currents in a 
Hodgkin-Huxley type neuron model, respectively (Ladenbauer 
et al., 2012). Furthermore, fitting the aEIF model parameters to 
a detailed biophysical model using standard electro-physiological 
paradigms revealed a clear relationship between parameter a and 
the conductance for the M current as well as between parameter 
b and the AHP current (not shown). 

Our method is based on several assumptions which allow to 
derive the mean-field equations. The Poisson approximation of 
spike train statistics is justified by experimental findings (Tolhurst 



et al, 1983; McAdams and Maunsell, 1999) although spiking 
seems to be more regular in some cortical areas (Maimon and 
Assad, 2009). The sparse random connectivity implies vanishing 
noise correlations between neurons in the large network limit and 
an experimental study in primary visual cortex of awake monkeys 
has reported almost zero noise correlations (Ecker et al., 2010). 
However, there is an ongoing debate about the strength of corre- 
lations in experimental data (Cohen and Kohn, 2011). We have 
used an adiabatic approximation, which relies on separable time 
scales of adaptation current and membrane voltage. Although this 
assumption is violated for small values of t w , numerically solv- 
ing the unreduced FP system, Equations (18)-(24), showed that 
our results are robust regarding the violation of this assumption. 
The results we obtained by simulations of aEIF networks and 
the mean-field results show quantitative differences. However, 
the presented effects described using the mean-field model are 
validated qualitatively by the network simulations. 

We have shown that spike-triggered adaptation provides a 
mechanism to generate spike rate oscillations in a low frequency 
range (alpha band and lower) if recurrent excitation is suffi- 
ciently strong. Increased subthreshold adaptation on the other 
hand does not contribute to this mechanism but rather damp- 
ens such oscillations. The type of adaptation current therefore 
strongly determines rhythmic activity in excitation dominated 
networks. The importance of activity-driven adaptation for slow 
oscillations is consistent with results from simulations of detailed 
(thalamo-)cortical spiking neuron network models (Bazhenov 
et al, 2002; Compte et al, 2003; Destexhe, 2009), mean-field 
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FIGURE 6 | Effects of adaptation on resonance properties of recurrent 
networks. Gain (top) and phase shift (bottom) of the spike rate re for 
networks with dominant recurrent excitation [g= 1.05, (A)] and inhibition 
[g = 1.5, (B)] as a function of the input frequency /. The gain is defined as the 
quotient of the oscillation amplitude in f£ for the input with frequency f and 
the amplitude for the lowest frequency (f min = 0.5 Hz). Adaptation parameter 
values are a = b = 0 (black), a = 5 nS, b = 0.05 nA, t w = 100 ms (dark red), 



and a = 5 nS, b = 0.05 nA, t w = 400 ms (orange). Delay distributions are 
identically parameterized (p^ = pj) with r f = 1.25 ms and id = 1.5 ms. The 
Poisson-input rates rf xt , fj xt each consist of a baseline rate plus a sinusoidal 
component of small amplitude (1 /1000th of the baseline) with frequency f. 
The baseline of rf xt is chosen to yield a steady-state spike rate ff of 50 Hz 
with constant input rate. The baseline rate of rj xt is chosen as explained in 
the Methods section. 



studies based on networks of excitatory neurons under the 
assumption sparse (Gigante et al., 2007b) and all-to-all con- 
nectivity (Nesse et al, 2008), as well as phenomenological rate 
models (Latham et al., 2000). We have further shown that reduc- 
ing inhibitory synaptic strength leads to a reduction on oscillation 
frequency, which is in agreement with similar experimental find- 
ings ( Sanchez- Vives et al, 2010). 

The M and AHP K + currents, which mediate spike frequency 
adaptation in pyramidal neurons, are known to be deactivated by 
acetylcholine (McCormick, 1992), with the AHP current show- 
ing higher sensitivity. Since the adaptation parameter b is strongly 
related to AHP type adaptation, our results support the hypoth- 
esis that the cholinergically induced activating transition from 
slow- wave oscillations to asynchronous irregular states (Lee and 
Dan, 2012) is mediated (at least in part) by a reduction of 
spike-triggered adaptation (Destexhe, 2009). 

We have demonstrated that an increase of either type of adap- 
tation current leads to a reduction in the frequency of oscillations 
generated by a loop of recurrent excitation and inhibition. This 
shows that the dynamical properties of neurons in addition to 
coupling characteristics strongly affect the network frequency. 
Also the passive (integrative) membrane properties significantly 
influence such networks oscillations as has been described previ- 
ously (Geisler et al., 2005). Our additional finding of decreased 



frequencies for increased adaptation time constants is consistent 
with the results from a computational study on clustering effects 
of spike-triggered adaptation in gamma oscillations (Kilpatrick 
and Ermentrout, 2011). 

Low input frequencies have been shown to be suppressed in the 
output of single excitatory neurons with increased spike-triggered 
(Gigante et al, 2007a) or subthreshold adaptation (Richardson 
et al., 2003; Prescott and Sejnowski, 2008), which we confirmed 
using our aEIF-based mean-field model. Such a high pass prop- 
erty of single neurons has also been found using a more general 
model of adaptation (Benda and Herz, 2003). We have demon- 
strated that both adaptation currents cause spike rate resonance 
in excitation dominated recurrent networks. Inhibition domi- 
nated networks, on the other hand, exhibit resonance without 
adaptation and we have shown that increased adaptation of exci- 
tatory neurons strongly amplifies this resonance. A similar effect 
has been described for purely inhibitory networks (Richardson, 
2009). In addition, our results show that adaptation shifts the 
resonance frequency to lower values. 

In excitation dominated networks, adaptation further leads to 
phase advances for low input frequencies in addition to phase 
delays for higher frequencies as observed in previous studies 
on single excitatory neurons (Fuhrmann et al., 2002; Gigante 
et al, 2007a). These adaptation-induced phase advances enable 
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synchronization of periodic activity between distant neurons (and 
populations of neurons) in different areas of the brain if the 
strength of adaptation is controlled appropriately, e.g., through 
cholinergic neuromodulation. 

Here we have considered one adaptation current for each neu- 
ron of the excitatory population. To account for the multimodal 
distribution of adaptation timescales found experimentally (La 



Camera et al., 2006) our approach can be easily extended to 
include multiple adaptation currents. 
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